﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using PaintDotNet;
using PaintDotNet.Effects;

namespace LibFiltr
{
    public class AnisotropyFilter
    {
        double _kappa;
        double dd, cdd;

        public AnisotropyFilter(double kappa = 30.0)
        {
            _kappa = kappa;
            dd = Math.Sqrt(2.0);
            cdd = 1 / Math.Pow(dd, 2);
        }

        private enum Channel
        {
            Red, Green, Blue
        }

        public void Filter(Surface dst, Surface src, Rectangle selection, Rectangle rect)
        {
            byte[,] r = FilterChannel(Channel.Red, src, rect);
            byte[,] g = FilterChannel(Channel.Green, src, rect);
            byte[,] b = FilterChannel(Channel.Blue, src, rect);
            ColorBgra CurrentPixel;
            for (int y = rect.Top; y < rect.Bottom; y++)
                for (int x = rect.Left; x < rect.Right; x++)
                {
                    CurrentPixel = src[x, y];
                    int i = y - rect.Top, j = x - rect.Left;
                    CurrentPixel.R = r[i, j];
                    CurrentPixel.G = g[i, j];
                    CurrentPixel.B = b[i, j];
                    dst[x, y] = CurrentPixel;
                }
        }

        private byte[,] FilterChannel(Channel ch, Surface src, Rectangle rect)
        {
            int[,] nablaN = CalcNablaN(ch, src, rect);
            int [,] nablaS = CalcNablaS(ch, src, rect);
            int [,] nablaW = CalcNablaW(ch, src, rect);
            int [,] nablaE = CalcNablaE(ch, src, rect);
            int [,] nablaNE = CalcNablaNE(ch, src, rect);
            int [,] nablaSE = CalcNablaSE(ch, src, rect);
            int [,] nablaSW = CalcNablaSW(ch, src, rect);
            int[,] nablaNW = CalcNablaNW(ch, src, rect);
            Console.WriteLine("all nablas calcled");

            double[,] cN = CalcCoeff(nablaN);
            double[,] cS = CalcCoeff(nablaS);
            double[,] cW = CalcCoeff(nablaW);
            double[,] cE = CalcCoeff(nablaE);
            double[,] cNE = CalcCoeff(nablaNE);
            double[,] cSE = CalcCoeff(nablaSE);
            double[,] cSW = CalcCoeff(nablaSW);
            double[,] cNW = CalcCoeff(nablaNW);
            Console.WriteLine("all coeffs calcled");

            byte[,] res = new byte[rect.Bottom - rect.Top, rect.Right - rect.Left];
            for (int y = rect.Top; y < rect.Bottom; y++)
            {
                for (int x = rect.Left; x < rect.Right; x++)
                {
                    int cv = GetChannel(src, x, y, ch);
                    int i = y - rect.Top, j = x - rect.Left;
                    double dv = cN[i, j] * nablaN[i, j] + cS[i, j] * nablaS[i, j] +
                                cW[i, j] * nablaW[i, j] + cE[i, j] * nablaE[i, j] +
                                cdd * cNE[i, j] * nablaNE[i, j] + cdd * cSE[i, j] * nablaSE[i, j] +
                                cdd * cSW[i, j] * nablaSW[i, j] + cdd * cNW[i, j] * nablaNW[i, j];
                    res[i, j] = (byte)Math.Round(cv + 0.1429 * dv);
                    Console.WriteLine(i.ToString() + " " + j.ToString());
                }
            }
            Console.WriteLine("bye :)");
            return res;
        }

        #region Nablas
        int[,] CalcNablaN(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top; y < rect.Bottom - 1; y++)
                for (int x = rect.Left; x < rect.Right; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x, y + 1, ch) - GetChannel(src, x, y, ch);
            // special case: last row
            int sy = rect.Bottom - 1;
            for (int x = rect.Left; x < rect.Right; x++)
                res[sy - rect.Top, x - rect.Left] = -GetChannel(src, x, sy, ch);
            return res;
        }

        int[,] CalcNablaS(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top + 1; y < rect.Bottom; y++)
                for (int x = rect.Left; x < rect.Right; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x, y - 1, ch) - GetChannel(src, x, y, ch);
            // special case: first row
            int sy = rect.Top;
            for (int x = rect.Left; x < rect.Right; x++)
                res[sy - rect.Top, x - rect.Left] = -GetChannel(src, x, sy, ch);
            return res;
        }

        int[,] CalcNablaW(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top; y < rect.Bottom; y++)
                for (int x = rect.Left; x < rect.Right - 1; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x + 1, y, ch) - GetChannel(src, x, y, ch);
            // special case: last column
            int sx = rect.Right - 1;
            for (int y = rect.Top; y < rect.Bottom; y++)
                res[y - rect.Top, sx - rect.Left] = -GetChannel(src, sx, y, ch);
            return res;
        }

        int[,] CalcNablaE(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top; y < rect.Bottom; y++)
                for (int x = rect.Left + 1; x < rect.Right; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x - 1, y, ch) - GetChannel(src, x, y, ch);
            // special case: first column
            int sx = rect.Left;
            for (int y = rect.Top; y < rect.Bottom; y++)
                res[y - rect.Top, sx - rect.Left] = -GetChannel(src, sx, y, ch);
            return res;
        }

        int[,] CalcNablaNE(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top; y < rect.Bottom - 1; y++)
                for (int x = rect.Left + 1; x < rect.Right; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x - 1, y + 1, ch) - GetChannel(src, x, y, ch);
            // special case: first column
            int sx = rect.Left;
            for (int y = rect.Top; y < rect.Bottom; y++)
                res[y - rect.Top, sx - rect.Left] = -GetChannel(src, sx, y, ch);
            // special case: last row
            int sy = rect.Bottom - 1;
            for (int x = rect.Left; x < rect.Right; x++)
                res[sy - rect.Top, x - rect.Left] = -GetChannel(src, x, sy, ch);
            return res;
        }

        int[,] CalcNablaSE(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top + 1; y < rect.Bottom; y++)
                for (int x = rect.Left + 1; x < rect.Right; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x - 1, y - 1, ch) - GetChannel(src, x, y, ch);
            // special case: first column
            int sx = rect.Left;
            for (int y = rect.Top; y < rect.Bottom; y++)
                res[y - rect.Top, sx - rect.Left] = -GetChannel(src, sx, y, ch);
            // special case: first row
            int sy = rect.Top;
            for (int x = rect.Left; x < rect.Right; x++)
                res[sy - rect.Top, x - rect.Left] = -GetChannel(src, x, sy, ch);
            return res;
        }

        int[,] CalcNablaSW(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top + 1; y < rect.Bottom; y++)
                for (int x = rect.Left; x < rect.Right - 1; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x + 1, y - 1, ch) - GetChannel(src, x, y, ch);
            // special case: last column
            int sx = rect.Right - 1;
            for (int y = rect.Top; y < rect.Bottom; y++)
                res[y - rect.Top, sx - rect.Left] = -GetChannel(src, sx, y, ch);
            // special case: first row
            int ly = rect.Top;
            for (int x = rect.Left; x < rect.Right; x++)
                res[ly - rect.Top, x - rect.Left] = -GetChannel(src, x, ly, ch);
            return res;
        }

        int[,] CalcNablaNW(Channel ch, Surface src, Rectangle rect)
        {
            int[,] res = new int[(rect.Bottom - rect.Top), (rect.Right - rect.Left)];
            for (int y = rect.Top; y < rect.Bottom - 1; y++)
                for (int x = rect.Left; x < rect.Right - 1; x++)
                    res[y - rect.Top, x - rect.Left] = GetChannel(src, x + 1, y + 1, ch) - GetChannel(src, x, y, ch);
            // special case: last column
            int sx = rect.Right - 1;
            for (int y = rect.Top; y < rect.Bottom; y++)
                res[y - rect.Top, sx - rect.Left] = -GetChannel(src, sx, y, ch);
            // special case: last row
            int ly = rect.Bottom - 1;
            for (int x = rect.Left; x < rect.Right; x++)
                res[ly - rect.Top, x - rect.Left] = -GetChannel(src, x, ly, ch);
            return res;
        }
        #endregion

        double[,] CalcCoeff(int[,] nabla)
        {
            int w = nabla.GetLength(1), h = nabla.GetLength(0);
            double[,] res = new double[h, w];
            for (int i = 0; i < h; i++)
                for (int j = 0; j < w; j++)
                    res[i, j] = Math.Exp(-Math.Pow(nabla[i,j] / _kappa, 2));
            return res;            
        }

        byte GetChannel(Surface s, int x, int y, Channel ch)
        {
            switch (ch)
            {
                case Channel.Red: return s[x, y].R;
                case Channel.Green: return s[x, y].G;
                case Channel.Blue: return s[x, y].B;
                default: return 0xFF;  // Make Compiler Happy =)
            }
        }
    }
}
